El CoVid-19 es una nueva enfermedad respiratoria causada por un virus. Entre sus síntomas podemos encontrar la fiebre la tos seca, además de otros problemas respiratorios. Lo peligroso de este virus es su facilidad para el contagio, pero una buena higiene y medidas como el confinamiento podrán detener su propagación.
Los datos que se van a analizar en este documento pertenecen a un dataset creado por usuarios de Kaggle.com. La fecha de este análisis comienza el día 24 de Abril de 2020, utilizando la versión número 87 del dataset proporcionado en el enlace anterior. El objetivo de este análisis es totalmente didáctico.
Fecha de última actualización: 24/04/2020
Testeamos una primera carga directamente con Python:
import pandas as pd
data = pd.read_csv("covid_19_clean_complete.csv")
#data.head(10)
Otra forma de testear la carga de datos, esta vez con código de R:
pd <- import("pandas")
data <- pd$read_csv("covid_19_clean_complete.csv")
#head(data)
Pero la librería que vamos a utilizar, es aquella que nos permite formatear las tablas y sacarlas de manera mucho más bonita y organizada. Para ello solo tenemos que escribir lo siguiente:
pd <- import("pandas")
data <- pd$read_csv("covid_19_clean_complete.csv")
kable(head(data))
| Province/State | Country/Region | Lat | Long | Date | Confirmed | Deaths | Recovered |
|---|---|---|---|---|---|---|---|
| NaN | Afghanistan | 33.0000 | 65.0000 | 1/22/20 | 0 | 0 | 0 |
| NaN | Albania | 41.1533 | 20.1683 | 1/22/20 | 0 | 0 | 0 |
| NaN | Algeria | 28.0339 | 1.6596 | 1/22/20 | 0 | 0 | 0 |
| NaN | Andorra | 42.5063 | 1.5218 | 1/22/20 | 0 | 0 | 0 |
| NaN | Angola | -11.2027 | 17.8739 | 1/22/20 | 0 | 0 | 0 |
| NaN | Antigua and Barbuda | 17.0608 | -61.7964 | 1/22/20 | 0 | 0 | 0 |
#kable(head(data, 10))
Aunque sí que es cierto que esta sintaxis puede resultar, a priori, algo confusa. Es por ello que vamos a utilizar otra librería para hacer todo mucho más “legible”. Esta librería no es otra que ‘tidyverse’. Para ello, usamos código nativo de R:
data <- read.csv("covid_19_clean_complete.csv", stringsAsFactors = FALSE)
data %>% head(10) %>% kable()
| Province.State | Country.Region | Lat | Long | Date | Confirmed | Deaths | Recovered |
|---|---|---|---|---|---|---|---|
| Afghanistan | 33.0000 | 65.0000 | 1/22/20 | 0 | 0 | 0 | |
| Albania | 41.1533 | 20.1683 | 1/22/20 | 0 | 0 | 0 | |
| Algeria | 28.0339 | 1.6596 | 1/22/20 | 0 | 0 | 0 | |
| Andorra | 42.5063 | 1.5218 | 1/22/20 | 0 | 0 | 0 | |
| Angola | -11.2027 | 17.8739 | 1/22/20 | 0 | 0 | 0 | |
| Antigua and Barbuda | 17.0608 | -61.7964 | 1/22/20 | 0 | 0 | 0 | |
| Argentina | -38.4161 | -63.6167 | 1/22/20 | 0 | 0 | 0 | |
| Armenia | 40.0691 | 45.0382 | 1/22/20 | 0 | 0 | 0 | |
| Australian Capital Territory | Australia | -35.4735 | 149.0124 | 1/22/20 | 0 | 0 | 0 |
| New South Wales | Australia | -33.8688 | 151.2093 | 1/22/20 | 0 | 0 | 0 |
Una vez tenemos cargado el dataset, el siguiente paso lógico es intentar entender correctamente el mismo. Para ello se suelen hacer una serie de operaciones típicas, como cambiar el nombre a las columnas, organizarlas, etc. Veamos rápidamente cómo es la estructura de los datos:
str(data)
## 'data.frame': 24366 obs. of 8 variables:
## $ Province.State: chr "" "" "" "" ...
## $ Country.Region: chr "Afghanistan" "Albania" "Algeria" "Andorra" ...
## $ Lat : num 33 41.2 28 42.5 -11.2 ...
## $ Long : num 65 20.17 1.66 1.52 17.87 ...
## $ Date : chr "1/22/20" "1/22/20" "1/22/20" "1/22/20" ...
## $ Confirmed : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Deaths : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Recovered : int 0 0 0 0 0 0 0 0 0 0 ...
Aquí podemos observar que estamos trabajando con un DataFrame, que es la estructura típica para trabajar en Python y R, que consta de 24366 observaciones en 8 columnas. Además, nos muestra el tipo de cada columna, cosa que es realmente úitl para saber con qué estamos trabajando en cada momento. Lo que sí es cierto, es que quizá los nombres de las columnas no son los adecuados para hacer un buen informe. Primero porque no están en el idioma en el cuál está hecho el informe, y segundo por el formato que tiene. Por ello, vamos a cambiarlos de manera que quede más legible y estético.
colnames(data)
## [1] "Province.State" "Country.Region" "Lat" "Long"
## [5] "Date" "Confirmed" "Deaths" "Recovered"
colnames(data) = c("Provincia_Estado",
"Pais_Region",
"Latitud", #Norte (+) o Sur (-)
"Longitud", #Oeste (-) o Este (+)
"Fecha",
"Casos_Confirmados", #Número de casos confirmados
"Casos_Muertos", #Número de muertes
"Casos_Recuperados") #Número de Recuperados
data %>% head(10) %>% kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados |
|---|---|---|---|---|---|---|---|
| Afghanistan | 33.0000 | 65.0000 | 1/22/20 | 0 | 0 | 0 | |
| Albania | 41.1533 | 20.1683 | 1/22/20 | 0 | 0 | 0 | |
| Algeria | 28.0339 | 1.6596 | 1/22/20 | 0 | 0 | 0 | |
| Andorra | 42.5063 | 1.5218 | 1/22/20 | 0 | 0 | 0 | |
| Angola | -11.2027 | 17.8739 | 1/22/20 | 0 | 0 | 0 | |
| Antigua and Barbuda | 17.0608 | -61.7964 | 1/22/20 | 0 | 0 | 0 | |
| Argentina | -38.4161 | -63.6167 | 1/22/20 | 0 | 0 | 0 | |
| Armenia | 40.0691 | 45.0382 | 1/22/20 | 0 | 0 | 0 | |
| Australian Capital Territory | Australia | -35.4735 | 149.0124 | 1/22/20 | 0 | 0 | 0 |
| New South Wales | Australia | -33.8688 | 151.2093 | 1/22/20 | 0 | 0 | 0 |
El siguiente problema que nos podemos encontrar es que alguna de las columnas no tenga el tipo de dato adecuado. Por ejemplo, tenemos datos cuantitativos, cualtitativos y ordinales (aquellos datos cualtitativos con los que, por ejemplo, podemos “calcular la media”).
A modo de ejemplo, la primera columna “Provincia_Estado”, describe una cualidad; es decir, es una variable cualitativa, al igual que la siguiente columna. El primer problema nos lo encontramos con la columna “Fecha”, ya que debería ser una variable ordinal, y no cualitativa (desde un punto de vista subjetivo, dado que también podría ser un factor ordenado y no una variable ordinal). Como vemos, podemos decir que las fechas son unas variables ciertamente “especiales”. Es por esto que existen varios paquetes con los que trabajarlas. Convirtamos pues, la variable “Fecha”, “Provincia_Estado” y “Pais_Region”:
data$Provincia_Estado %<>% factor()
data$Pais_Region %<>% factor()
#data$Fecha %<>% as.Date(format="%m/%d/%y") #Prueba de formateo de fecha con R.
data$Fecha %<>% mdy() #Uso de Lubridate para la conversión
str(data)
## 'data.frame': 24366 obs. of 8 variables:
## $ Provincia_Estado : Factor w/ 81 levels "","Alberta","Anguilla",..: 1 1 1 1 1 1 1 1 6 49 ...
## $ Pais_Region : Factor w/ 185 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 9 ...
## $ Latitud : num 33 41.2 28 42.5 -11.2 ...
## $ Longitud : num 65 20.17 1.66 1.52 17.87 ...
## $ Fecha : Date, format: "2020-01-22" "2020-01-22" ...
## $ Casos_Confirmados: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Casos_Muertos : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Casos_Recuperados: int 0 0 0 0 0 0 0 0 0 0 ...
Como vemos ahora tenemos todo en un formato correcto. El uso de lubridate trae muchas ventajas a nuestro entorno de trabajo. Una de ellas, por ejemplo, son los cálculos rápidos de variables que en un principio, son categóricas ordenadas:
dmy("06/04/20") -> d2
dmy("21/01/20") -> d1
d2-d1
## Time difference of 76 days
Otro problema con el que nos encontramos a la hora de hacer un data wrangling, es que es posible que exitan datos que no tengan demasiado sentido, que sean anómalos. Una comprobación simple que podríamos hacer es, sabiendo que:
\[ Casos \hspace{0.1cm} confirmados \hspace{0.1cm} = Casos \hspace{0.1cm} muertos \hspace{0.1cm}+ \hspace{0.1cm} recuperados \hspace{0.1cm}+\hspace{0.1cm}enfermos \] Calcular el número de Enfermos, por si en cualquier momento necesitarámos saber este dato. Para ello, podemos usar la función mutate de tidyverse:
data %<>%
mutate(Casos_Enfermos = Casos_Confirmados - Casos_Muertos - Casos_Recuperados)
data %>%
filter(Casos_Confirmados > 10000) %>%
head() %>%
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|---|---|---|---|
| Hubei | China | 30.9756 | 112.2707 | 2020-02-02 | 11177 | 350 | 295 | 10532 |
| Hubei | China | 30.9756 | 112.2707 | 2020-02-03 | 13522 | 414 | 386 | 12722 |
| Hubei | China | 30.9756 | 112.2707 | 2020-02-04 | 16678 | 479 | 522 | 15677 |
| Hubei | China | 30.9756 | 112.2707 | 2020-02-05 | 19665 | 549 | 633 | 18483 |
| Hubei | China | 30.9756 | 112.2707 | 2020-02-06 | 22112 | 618 | 817 | 20677 |
| Hubei | China | 30.9756 | 112.2707 | 2020-02-07 | 24953 | 699 | 1115 | 23139 |
Con esta columna, tenemos la oportunidad de comprobar si existe algún dato anómalo, ya que siguiendo la ecuación anterior, jamás la suma de casos muertos, recuperados y enfermos debe superar al número de casos confirmados. Para esto podemos hacer un filtrado simple.
data %>%
filter(Casos_Enfermos < 0) %>%
arrange(Provincia_Estado, Fecha) %>% #Para ordenar por Provincia_Estado y Fecha
head() %>%
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|---|---|---|---|
| Diamond Princess | Canada | 0 | 0 | 2020-03-22 | 0 | 1 | 0 | -1 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-23 | 0 | 1 | 0 | -1 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-24 | 0 | 1 | 0 | -1 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-25 | 0 | 1 | 0 | -1 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-26 | 0 | 1 | 0 | -1 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-27 | 0 | 1 | 0 | -1 |
Oh sorpresa, sí que hay unos cuantos. Alguien está mintiendo con sus datos o se ha equivocado, ya que el de Canadá tiene latitud y longitud 0, cosa imposible. Por seguir comprobando, vamos a ver qué pasó en Hainan.
data %>%
filter(Provincia_Estado == "Hainan") %>%
head() %>% #Borrar esta línea para visualización completa
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|---|---|---|---|
| Hainan | China | 19.1959 | 109.7453 | 2020-01-22 | 4 | 0 | 0 | 4 |
| Hainan | China | 19.1959 | 109.7453 | 2020-01-23 | 5 | 0 | 0 | 5 |
| Hainan | China | 19.1959 | 109.7453 | 2020-01-24 | 8 | 0 | 0 | 8 |
| Hainan | China | 19.1959 | 109.7453 | 2020-01-25 | 19 | 0 | 0 | 19 |
| Hainan | China | 19.1959 | 109.7453 | 2020-01-26 | 22 | 0 | 0 | 22 |
| Hainan | China | 19.1959 | 109.7453 | 2020-01-27 | 33 | 1 | 0 | 32 |
Vemos pues, que el posible error para esta provincia es que es posible que dejaran de contar. Ya que vemos que ese 168 se mantiene constante, hasta el 1 de Abril de 2020, donde lo corrigen sumando el -6 a la columna de casos recuperados. Podemos tratar de corregir estos datos de la siguiente forma:
data %>%
filter(Provincia_Estado == "Hainan", Casos_Enfermos < 0) %>%
mutate(Casos_Recuperados = Casos_Recuperados + Casos_Enfermos,
Casos_Enfermos = 0) %>%
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|---|---|---|---|
| Hainan | China | 19.1959 | 109.7453 | 2020-03-24 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-25 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-26 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-27 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-28 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-29 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-30 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-03-31 | 168 | 6 | 162 | 0 |
| Hainan | China | 19.1959 | 109.7453 | 2020-04-01 | 168 | 6 | 162 | 0 |
En cambio, los datos de Canadá parecen ser una prueba que se olvidaron quitar. Para que no afecte, cambiemos este número de la misma forma que hemos hecho en China, que es algo que no afectará a nuestro resultado.
data %>%
filter(Provincia_Estado == "Diamond Princess", Casos_Enfermos < 0) %>%
mutate(Casos_Recuperados = Casos_Recuperados + Casos_Enfermos,
Casos_Enfermos = 0) %>%
head() %>%
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|---|---|---|---|
| Diamond Princess | Canada | 0 | 0 | 2020-03-22 | 0 | 1 | -1 | 0 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-23 | 0 | 1 | -1 | 0 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-24 | 0 | 1 | -1 | 0 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-25 | 0 | 1 | -1 | 0 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-26 | 0 | 1 | -1 | 0 |
| Diamond Princess | Canada | 0 | 0 | 2020-03-27 | 0 | 1 | -1 | 0 |
Si quisiéramos por ejemplo, analizar qué ocurre en Europa, no tenemos una forma rápida y eficiente de saber, con los datos que tenemos en nuestro dataset, qué países pertenecen o no a Europa. El truco que aquí utilizaré es encuadrar en un mapamundi Europa y apuntar los valores máximos y mínimos de latitud y longitud de manera que podamos hacer un filtrado.
#datos_europa = data[data$Latitud >38 & data$Longitud > -25 & data$Longitud < 30,]
datos_europa = data %>%
filter(Latitud > 38, between(Longitud, -25, 30))
nrow(datos_europa)
## [1] 4185
table(datos_europa$Pais_Region) %>%
as.data.frame() %>%
filter(Freq > 0) %>%
kable()
| Var1 | Freq |
|---|---|
| Albania | 93 |
| Andorra | 93 |
| Austria | 93 |
| Belarus | 93 |
| Belgium | 93 |
| Bosnia and Herzegovina | 93 |
| Bulgaria | 93 |
| Croatia | 93 |
| Czechia | 93 |
| Denmark | 186 |
| Estonia | 93 |
| Finland | 93 |
| France | 93 |
| Germany | 93 |
| Greece | 93 |
| Holy See | 93 |
| Hungary | 93 |
| Iceland | 93 |
| Ireland | 93 |
| Italy | 93 |
| Kosovo | 93 |
| Latvia | 93 |
| Liechtenstein | 93 |
| Lithuania | 93 |
| Luxembourg | 93 |
| Moldova | 93 |
| Monaco | 93 |
| Montenegro | 93 |
| Netherlands | 93 |
| North Macedonia | 93 |
| Norway | 93 |
| Poland | 93 |
| Portugal | 93 |
| Romania | 93 |
| San Marino | 93 |
| Serbia | 93 |
| Slovakia | 93 |
| Slovenia | 93 |
| Spain | 93 |
| Sweden | 93 |
| Switzerland | 93 |
| United Kingdom | 279 |
En este caso parece que hemos hecho un buen filtrado. Aunque no siempre se suele hacer así, si no mediante un filtrado por bolas. Esto es, en vez de coger un rectángulo mediante latitud y longitud, cogemos un círculo que esté centrado en el sitio que queremos, ya sea Madrid, Potsdam, o Murcia. Básicamente mediante la distancia euclídea.
\[d(x,y) = \sqrt{(x_{Lat}-y_{Lat})^2+(x_{Long}-y_{Long})^2}\] Por lo tanto, podemos definir una función distancia del siguiente modo:
distancia_grados = function(x,y){
sqrt((x[1]-y[1])^2+(x[2]-y[2])^2)
}
Ahora podemos calcular, por ejemplo, la distancia a Madrid.
distancia_grados_madrid = function(x){
madrid = c(40.4378693,-3.8199644) #Coordenadas obtenidas de Maps.
distancia_grados(x, madrid)
}
Entonces ahora, podemos coger todo el dataset de, por ejemplo, Europa, y calcular la distancia a Madrid.
dist_madrid = apply(cbind(datos_europa$Latitud, datos_europa$Longitud),
MARGIN = 1, #Para Filas
FUN = distancia_grados_madrid) #Aplicamos la función
datos_europa %<>%
mutate(dist_madrid = dist_madrid)
datos_europa %>%
filter(between(Fecha, dmy("02/03/20"), dmy("10/03/20")), #Por poner unas fechas
dist_madrid < 10) %>% #Una bola de radio 10 grados
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos | dist_madrid |
|---|---|---|---|---|---|---|---|---|---|
| Andorra | 42.5063 | 1.5218 | 2020-03-02 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-02 | 191 | 3 | 12 | 176 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-02 | 2 | 0 | 0 | 2 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-02 | 120 | 0 | 2 | 118 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-02 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-03 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-03 | 204 | 4 | 12 | 188 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-03 | 2 | 0 | 0 | 2 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-03 | 165 | 1 | 2 | 162 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-03 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-04 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-04 | 285 | 4 | 12 | 269 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-04 | 5 | 0 | 0 | 5 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-04 | 222 | 2 | 2 | 218 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-04 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-05 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-05 | 377 | 6 | 12 | 359 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-05 | 8 | 0 | 0 | 8 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-05 | 259 | 3 | 2 | 254 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-05 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-06 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-06 | 653 | 9 | 12 | 632 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-06 | 13 | 0 | 0 | 13 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-06 | 400 | 5 | 2 | 393 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-06 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-07 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-07 | 949 | 11 | 12 | 926 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-07 | 20 | 0 | 0 | 20 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-07 | 500 | 10 | 30 | 460 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-07 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-08 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-08 | 1126 | 19 | 12 | 1095 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-08 | 30 | 0 | 0 | 30 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-08 | 673 | 17 | 30 | 626 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-08 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-09 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-09 | 1209 | 19 | 12 | 1178 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-09 | 30 | 0 | 0 | 30 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-09 | 1073 | 28 | 32 | 1013 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-09 | 0 | 0 | 0 | 0 | 9.0522218 |
| Andorra | 42.5063 | 1.5218 | 2020-03-10 | 1 | 0 | 0 | 1 | 5.7282504 | |
| France | 46.2276 | 2.2137 | 2020-03-10 | 1784 | 33 | 12 | 1739 | 8.3621820 | |
| Portugal | 39.3999 | -8.2245 | 2020-03-10 | 41 | 0 | 0 | 41 | 4.5251866 | |
| Spain | 40.0000 | -4.0000 | 2020-03-10 | 1695 | 35 | 32 | 1628 | 0.4734367 | |
| Channel Islands | United Kingdom | 49.3723 | -2.3644 | 2020-03-10 | 1 | 0 | 0 | 1 | 9.0522218 |
De esta forma podemos analizar de manera fácil durante las fechas propuestas, cuántos casos teníamos en torno a Madrid. Aunque hay que tener cierto cuidado, dado que estamos hablando en términos de países, no solo de Madrid. Siempre hay que saber qué se está analizando.
Pero, ¿qué es un análisis geográfico sin su correspondiente representación? Para ello vamos a utilizar las librerías ggmaps y rnaturalearth.
world <- ne_countries(scale = "medium", returnclass = "sf") #Definimos un mapa
ggplot(data = world) + #Especificamos las capas de dibujo
geom_sf() + #Aquí pdoríamos añadir color="black", fill="green" por ejemplo
xlab("Longitud") + ylab("Latitud") +
ggtitle("Mapa del Mundo", subtitle = "CoVid-19")
También podemos retocar la estética del mapa. Podemos hacer algo como:
ggplot(data = world) +
geom_sf(color = "black", aes(fill = mapcolor13)) +
xlab("Longitud") + ylab("Latitud") +
ggtitle("Mapa del Mundo", subtitle = "CoVid-19")
Es decir, podemos ver que este paquete ya tiene un montón de información cuando creamos world. Para poder sacar partido a esto, tenemos que cruzar los datos de nuestro dataset con la información que nos ha creado este paquete. Esto puede ser un jaleo terrible ya que los nombres pueden o no coincidir, etc.
world <- ne_countries(scale = "medium", returnclass = "sf")
world %>%
inner_join(data, by = c("name" = "Pais_Region")) %>%
filter(Fecha == dmy("15/03/2020")) %>%
ggplot() +
geom_sf(color = "black", aes(fill = Casos_Confirmados)) +
xlab("Longitud") + ylab("Latitud") +
ggtitle("Mapa del Mundo", subtitle = "CoVid-19")
## Warning: Column `name`/`Pais_Region` joining character vector and factor,
## coercing into character vector
Como vemos, fallan países como los Estados Unidos, o el centro de África. Por lo tanto, podemos ver que esto no es que sea un método muy efectivo ya que a lo mejor los datos no están pensados precisamente para ser trabajados con este paquete. Pero sí que es posible que consigamos algo más minimalista, representando la concentración mediante por ejemplo, unas bolas. Para ello, podemos hacer lo siguiente:
data %>%
filter(Fecha == dmy("30/03/2020")) %>%
ggplot(aes(Longitud, Latitud)) +
geom_point(aes(size = Casos_Confirmados, colour = Casos_Muertos)) +
coord_fixed() +
theme(legend.position = "bottom")
Estos mapas también se pueden convertir en interactivos mediante el uso de Plotly. Pero en realidad los mapas, como ya hemos dicho, no aportan demasiado. Podemos hacer para tener un estudio completo, un top que nos permita visualizar por ejemplo cuál es la proporción de muertos, o el número de infectados.
data %>%
filter(Fecha == ymd("2020-04-05"),
Casos_Confirmados > 1000) %>%
mutate(Prop_Muertos = Casos_Muertos / Casos_Confirmados,
Ranking = dense_rank(desc(Prop_Muertos))) %>%
arrange(Ranking) %>%
head() %>%
kable()
| Provincia_Estado | Pais_Region | Latitud | Longitud | Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos | Prop_Muertos | Ranking |
|---|---|---|---|---|---|---|---|---|---|---|
| Italy | 43.0000 | 12.0000 | 2020-04-05 | 128948 | 15887 | 21815 | 91246 | 0.1232047 | 1 | |
| Algeria | 28.0339 | 1.6596 | 2020-04-05 | 1320 | 152 | 90 | 1078 | 0.1151515 | 2 | |
| France | 46.2276 | 2.2137 | 2020-04-05 | 70478 | 8078 | 16183 | 46217 | 0.1146173 | 3 | |
| United Kingdom | 55.3781 | -3.4360 | 2020-04-05 | 47806 | 4934 | 135 | 42737 | 0.1032088 | 4 | |
| Netherlands | 52.1326 | 5.2913 | 2020-04-05 | 17851 | 1766 | 250 | 15835 | 0.0989300 | 5 | |
| Spain | 40.0000 | -4.0000 | 2020-04-05 | 131646 | 12641 | 38080 | 80925 | 0.0960227 | 6 |
Para este tipo de datos podemos hacer un tratamiento con mosaic plots, algo muy interesante sobretodo para saber manejar en un futuro.
data$lat_class = cut(data$Latitud, breaks = seq(from = -90, to = 90, by = 10))
data$long_class = cut(data$Longitud, breaks = seq(from = -180, to = 180, by = 10))
tt = table(data$lat_class, data$long_class)
tt = tt[nrow(tt):1, ]
mosaicplot(t(tt), shade = TRUE, xlab = NULL, ylab = NULL)
Aunque las leyendas en este gráfico no le hacen mucho honor, porque son demasiados datos y no se lee, al fin y al cabo lo que estamos viendo aquí es un mapa del mundo, marcado por longitudes y latitudes. Donde Europa estaría entorno al centro, América a la izquierda, etc. Dónde los cuadrados son proporcionales al número de observaciones que tenemos.
El análiss temporal es muy utilizado en temas de enconomía, y es algo más complejo de lo que puedo mostrar aquí. Da para muchas horas de investigación. Así que aquí haremos algo sencillo para los datos que tenemos. A la hora de tartar de predecir, no siempre los datos van a ser lo aceptables que se podría desear. Es por esto que ahora voy a realizar un análisis meramente descriptivo.
datos_por_fecha = aggregate( #Para agregar los valores por fecha
cbind(Casos_Confirmados, Casos_Muertos, Casos_Recuperados, Casos_Enfermos) ~ Fecha,
data = data, #DataSet completo
FUN = sum #Función de agregación a utilizar
)
datos_por_fecha %>% head() %>% kable()
| Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|
| 2020-01-22 | 555 | 17 | 28 | 510 |
| 2020-01-23 | 654 | 18 | 30 | 606 |
| 2020-01-24 | 941 | 26 | 35 | 880 |
| 2020-01-25 | 1434 | 42 | 38 | 1354 |
| 2020-01-26 | 2118 | 56 | 51 | 2011 |
| 2020-01-27 | 2927 | 82 | 58 | 2787 |
Podemos hacer un par de representaciones para ver esto mejor:
#barplot(Casos_Confirmados ~ Fecha, data = datos_por_fecha) #Comprobación
fig <- plot_ly( #Figura Principal
x = datos_por_fecha$Fecha,
y = datos_por_fecha$Casos_Muertos,
name = "Casos Muertos",
type = "bar")
fig <- fig %>% add_trace( #Trazas adicionales
y = datos_por_fecha$Casos_Enfermos,
name = 'Casos Enfermos')
fig <- fig %>% add_trace(
y = datos_por_fecha$Casos_Recuperados,
name = 'Casos Recuperados')
#Damos formato a la figura
fig <- fig %>% layout(title = 'Análisis temporal' ,yaxis = list(title = 'Número de Casos'), xaxis = list(title = 'Fecha'), barmode = 'stack')
fig
fig <- plot_ly(
datos_por_fecha,
x = ~datos_por_fecha$Fecha,
y = ~datos_por_fecha$Casos_Confirmados,
name = 'Casos Confirmdos',
type = 'scatter',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~datos_por_fecha$Casos_Muertos,
name = 'Casos Muertos',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~datos_por_fecha$Casos_Recuperados,
name = 'Casos Recuperados',
color = I("green"),
mode = 'lines')
#Damos formato a la figura
fig <- fig %>% layout(title = 'Análisis temporal' ,yaxis = list(title = 'Número de Casos'), xaxis = list(title = 'Fecha'), barmode = 'stack')
fig
Aquí ya se puede ver como al priuncipio solo documentaban casos en China, pero a medida que pasa el tiempo y hasta que se declara una pandemia mundial, el número de casos crece con esa tendencia exponencial de la que tanto se habla.
Estos gráficos están bien, pero pueden llevar a cierta confusíon, ya que al estar las barras apiladas puede que se muestre una tendencia mucho más pronunciada de lo que es realmente. Para realizar un buen análisis podemos utilizar la librería xts, la cuál es muy útil para el análisis temporal.
datos_por_fecha_ts = xts(x = datos_por_fecha[, 2:5],
order.by = datos_por_fecha$Fecha)
dygraph(datos_por_fecha_ts) %>%
dyCSS("legend.css") %>%
dyOptions(labelsUTC = TRUE, labelsKMB = TRUE,
fillGraph = TRUE, fillAlpha = 0.1,
drawGrid = FALSE) %>%
dyRangeSelector() %>%
dyCrosshair(direction = "both") %>%
dyHighlight(highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE) %>%
dyRoller(rollPeriod = 2)
Podemos hacer una foto de la situación de lo que ocurría en España a lo largo del tiempo mediante un filtrado. Además, podemos utilizar la librería plotly para hacer alguna representación. Por ejemplo:
datos_españa1 = data %>%
filter(Pais_Region == "Spain", Fecha > dmy("01/03/20")) %>%
select(Fecha, starts_with("Casos_")) #Filtrado por columnas
datos_españa1 %>% head() %>% kable()
| Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos |
|---|---|---|---|---|
| 2020-03-02 | 120 | 0 | 2 | 118 |
| 2020-03-03 | 165 | 1 | 2 | 162 |
| 2020-03-04 | 222 | 2 | 2 | 218 |
| 2020-03-05 | 259 | 3 | 2 | 254 |
| 2020-03-06 | 400 | 5 | 2 | 393 |
| 2020-03-07 | 500 | 10 | 30 | 460 |
#Primera Figura del SubPlot
fig1 <- plot_ly(width = 1000, height = 400,
x = ~datos_españa1$Fecha,
y = ~datos_españa1$Casos_Confirmados,
type = 'scatter', mode = 'lines',
name = 'Casos Confirmados', fill = 'tozeroy')
fig1 <- fig1 %>% add_trace(x = ~datos_españa1$Fecha,
y = ~datos_españa1$Casos_Enfermos,
name = 'Casos Enfermos',
opacity = 0.4)
fig1 <- fig1 %>% add_trace(x = ~datos_españa1$Fecha,
y = ~datos_españa1$Casos_Recuperados,
name = 'Casos Recuperados',
opacity = 0.8)
fig1 <- fig1 %>% add_trace(x = ~datos_españa1$Fecha,
y = ~datos_españa1$Casos_Muertos,
name = 'Casos Muertos',
opacity = 0.8)
fig1 <- fig1 %>% layout(
xaxis = list(title = 'Fecha'),
yaxis = list(title = 'Número de Casos'))
#Segunda Figura del SubPlot
fig2 <- plot_ly(width = 1000, height = 400,
x = datos_españa1$Fecha,
y = datos_españa1$Casos_Enfermos,
name = "Casos Enfermos",
color = I("orange"),
type = "bar")
fig2 <- fig2 %>% add_trace(
y = datos_españa1$Casos_Recuperados,
name = 'Casos Recuperados',
color = I("green"))
fig2 <- fig2 %>% add_trace(
y = datos_españa1$Casos_Muertos,
color = I("red"),
name = 'Casos Muertos')
fig2 <- fig2 %>% layout(yaxis = list(title = 'Número de Casos'),
xaxis = list(title = 'Fecha'),
barmode = 'stack')
#Creación del SubPlot
fig <- subplot(fig1, fig2, margin=0.05, shareY = TRUE)
fig <- fig %>% layout(title = "Número de Casos totales en España")
# Mostramos el resultado
fig
Es importante destacar que ahora mismo estamos estudiando valores acumulados; es decir, si vemos por ejemplo, los casos confirmados en España:
datos_españa2 = data %>%
filter(Pais_Region == "Spain") %>%
select(Casos_Confirmados)
datos_españa2 %>% head(15) %>% kable()
| Casos_Confirmados |
|---|
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 0 |
| 1 |
| 1 |
| 1 |
| 1 |
| 1 |
Vemos que hay un 1 que se repite, esto siempre es el mismo individuo, entonces dependiendo del análisis, podríamos estar analizando a la misma persona repetidas veces. Para ello, podemos hacer una desagregación. Para ello, podemos utilizar la función Lag, la cual hará que vayamos un día retrasados en el tiempo:
#lag(datos_españa2$Casos_Confirmados, n = 1)
Como vemos, ahora tenemos un NA en el primer día, esto es porque como he dicho, vamos un día retrasados. Lo mismo podemos hacer pero adelantando un día todo, para ello podemos usar Lead.
#lead(datos_españa2$Casos_Confirmados, n = 1)
Entonces ahora mediante una simple resta podemos tener el número de casos diarios. Por ejemoplo:
#datos_españa2$Casos_Confirmados - lag(datos_españa2$Casos_Confirmados, n = 1)
datos_españa = data %>%
filter(Pais_Region == "Spain") %>%
select(Fecha, starts_with("Casos_"))
datos_españa %<>%
mutate(Nuevos_Casos_Confirmados = Casos_Confirmados - lag(Casos_Confirmados, n=1),
Nuevos_Casos_Muertos = Casos_Muertos - lag(Casos_Muertos, n=1),
Nuevos_Casos_Recuperados = Casos_Recuperados - lag(Casos_Recuperados, n=1))
#datos_españa %>% head(15) %>% kable()
Evidentemente, el primer día no aparece, pero tampoco nos interesa. Ahora mismo lo que nos interesa de verdad es que tenemos exactamente el número de nuevos confirmados. Podemos hacer una representación simple para ver ese comportamiento, ahora, más caótico.
fig <- plot_ly(
datos_españa,
x = ~datos_españa$Fecha,
y = ~datos_españa$Nuevos_Casos_Confirmados,
name = 'Nuevos Casos Confirmdos',
type = 'scatter',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~datos_españa$Nuevos_Casos_Muertos,
name = 'Nuevos Casos Muertos',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~datos_españa$Nuevos_Casos_Recuperados,
name = 'Nuevos Casos Recuperados',
mode = 'lines')
#Damos formato a la figura
fig <- fig %>% layout(title = 'Análisis temporal' ,yaxis = list(title = 'Número de Nuevos Casos'), xaxis = list(title = 'Fecha'))
fig
También podemos calcular la Tasa de Variación Media, esta se define como (Casos en el día t - Casos en el día t-1)/Casos en el día t. De manera que aplicando esto para los casos confirmados:
datos_españa %<>%
mutate(TVM = (Casos_Confirmados - lag(Casos_Confirmados, n=1))/Casos_Confirmados)
#datos_españa %>% head(20) %>% kable()
fig <- plot_ly(
datos_españa,
x = ~datos_españa$Fecha,
y = ~datos_españa$TVM,
name = 'Tasa de Variación Media',
type = 'scatter',
mode = 'lines')
fig <- fig %>% layout(title = 'Tasa de Variación Media' ,yaxis = list(title = 'TVM (t)'), xaxis = list(title = 'Fecha'))
fig
Si queremos comparar cómo ha avanzado (temporalmente) el virus en varios países nos encontramos con un problema muy básico, y es que mientras que en un país ha empezado, por ejemplo, hace una semana, en otro puede llevar ya un mes o dos. Es por esto que la comparación no se puede efectuar de manera directa. Esa comparación implica tener que volver atrás en el tiempo a donde empezó todo. A esto es a lo que se le llama un análisis por cohortes. El cohorte no es más que el segmento. Entonces, en este caso, cada país volverá a la fecha en la que se originó el virus en el propio sitio, de manera que podamos escalar todas las fechas. Es decir, queremos hallar el día 0 de cada país, el último día que el mismo fue ‘normal’ en términos de coronavirus.
primer_contagio = data %>%
group_by(Pais_Region) %>%
filter(Casos_Confirmados > 0) %>%
summarise(Primer_Contagio = min(Fecha)-1) #Restamos 1 para obtener el día 0.
primer_contagio %>% head() %>% kable()
| Pais_Region | Primer_Contagio |
|---|---|
| Afghanistan | 2020-02-23 |
| Albania | 2020-03-08 |
| Algeria | 2020-02-24 |
| Andorra | 2020-03-01 |
| Angola | 2020-03-19 |
| Antigua and Barbuda | 2020-03-12 |
De esta manera tenemos para cada país la fecha del día 0, lo que buscábamos. Ahora podemos crear un dataset que nos muestre todos los días que han pasado desde ese primer contagio y que me resuma la información que luego utilizaré para representar en un diagrama de cohortes.
data_first = data %>%
inner_join(primer_contagio, by = "Pais_Region") %>%
mutate(Dias_Desde_PC = as.numeric(Fecha - Primer_Contagio)) %>%
filter(Dias_Desde_PC >= 0) %>%
group_by(Dias_Desde_PC, Pais_Region) %>% #Para poder hacer una mejor representación
summarise(Casos_Confirmados = sum(Casos_Confirmados),
Casos_Muertos = sum(Casos_Muertos),
Casos_Recuperados = sum(Casos_Recuperados),
Casos_Enfermos = sum(Casos_Enfermos))
#Representacion
data_first %>%
filter(Pais_Region %in% c("Spain", "Italy", "China", "US", "Germany")) %>%
ggplot(aes(x = Dias_Desde_PC, y = Casos_Confirmados)) +
geom_line(aes(col = Pais_Region)) +
xlab("Días desde el Primer Contagio") +
ylab("Número de Casos Confirmados") +
ggtitle("Análisis por Cohortes") +
theme(legend.position = "top", legend.title = element_blank()) -> g
ggplotly(g)
Nos faltan, por otra parte, los Casos Recuperados, Enfermos, y Muertos. Para ello, podemos representar los 4 gráficos.
data_first %>%
filter(Pais_Region %in% c("Spain", "Italy", "China", "US", "Germany")) %>%
ggplot(aes(x = Dias_Desde_PC, y = Casos_Confirmados)) +
geom_line(aes(col = Pais_Region)) +
xlab("Días desde el Primer Contagio") +
ylab("nº de Casos Confirmados") -> g_con
data_first %>%
filter(Pais_Region %in% c("Spain", "Italy", "China", "US", "Germany")) %>%
ggplot(aes(x = Dias_Desde_PC, y = Casos_Enfermos)) +
geom_line(aes(col = Pais_Region)) +
xlab("Días desde el Primer Contagio") +
ylab("nº de Casos Enfermos") -> g_en
data_first %>%
filter(Pais_Region %in% c("Spain", "Italy", "China", "US", "Germany")) %>%
ggplot(aes(x = Dias_Desde_PC, y = Casos_Recuperados)) +
geom_line(aes(col = Pais_Region)) +
xlab("Días desde el Primer Contagio") +
ylab("nº de Casos Recuperados") -> g_rec
data_first %>%
filter(Pais_Region %in% c("Spain", "Italy", "China", "US", "Germany")) %>%
ggplot(aes(x = Dias_Desde_PC, y = Casos_Muertos)) +
geom_line(aes(col = Pais_Region)) +
xlab("Días desde el Primer Contagio") +
ylab("nº de Casos Muertos") -> g_mu
figure <- ggarrange(g_con, g_en, g_rec, g_mu,
ncol = 2, nrow = 2,
common.legend = TRUE, legend = "bottom")
annotate_figure(figure, top = text_grob("Analisis por Cohortes"))
Ahora podemos visualizar de mejor manera nuestro análisis por cohortes.
Predecir como va a avanzar la epidemia es más complicado de lo que se puede decir de manera directa. Esto es debido a que los modelos pandémicos no son lineales, o se rigen por ecuaciones diferenciales. En esencia vamos a tener la variable independiente (x), que será el número de días desde el origen de la pandemia, y la variable dependiente, que será aquello que queremos predecir (y).
La idea es construir un modelo de regresión que sea capaz de predecir esa variable dependiente. Centrémonos por ejemplo, en España.
datos_españa$Dias = as.numeric(datos_españa$Fecha - dmy("22/01/2020"))
datos_españa %>% head() %>% kable()
| Fecha | Casos_Confirmados | Casos_Muertos | Casos_Recuperados | Casos_Enfermos | Nuevos_Casos_Confirmados | Nuevos_Casos_Muertos | Nuevos_Casos_Recuperados | TVM | Dias |
|---|---|---|---|---|---|---|---|---|---|
| 2020-01-22 | 0 | 0 | 0 | 0 | NA | NA | NA | NA | 0 |
| 2020-01-23 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | NaN | 1 |
| 2020-01-24 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | NaN | 2 |
| 2020-01-25 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | NaN | 3 |
| 2020-01-26 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | NaN | 4 |
| 2020-01-27 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | NaN | 5 |
Hay muchos tipos de modelos de predicción, empecemos por los más sencillos.
En R existe la función lm, que nos construye un modelo lineal de manera muy fácil:
mod1 <- lm(Casos_Confirmados ~ Dias, datos_españa)
summary(mod1)
##
## Call:
## lm(formula = Casos_Confirmados ~ Dias, data = datos_españa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57050 -35004 1801 34308 61433
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -54465.3 7811.2 -6.973 4.82e-10 ***
## Dias 2239.7 146.7 15.272 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37970 on 91 degrees of freedom
## Multiple R-squared: 0.7193, Adjusted R-squared: 0.7162
## F-statistic: 233.2 on 1 and 91 DF, p-value: < 2.2e-16
La ordenada en el origen tiene un p-valor muy muy pequeño, por lo tanto tiene mucha significatividad en el modelo. Además nos da el valor de \(R^{2}\). Que no está nada mal en realidad. Entonces, nuestro modelo se puede resumir en:
\[Casos\ Confirmados = 2239.744457 Dias + -5.4465288\times 10^{4}\]
Pero, ¿es nuestro modelo fiable? Hagamos una primera representación:
plot(datos_españa$Dias, datos_españa$Casos_Confirmados)
abline(mod1, col = "red")
Como vemos, parece el ajuste hecho por un mono. Podemos pintar también la distribución de los errores:
plot(mod1$residuals ~ mod1$fitted.values, xlab = "Valores Ajustados",
ylab = "Residuos del Modelo")
Aquí ya podemos observar dos patrones claramente diferenciados. Obviamente, si este modelo fuera bueno, no se observaría quí ningún patrón. Como preveíamos, este modelo tiene pinta de ser pésimo para este caso. Otra cosa que debería cumplirse, es que el modelo tuviera la misma varianza de los errores. Para comprobarlo hagamos lo siguiente:
residuos = mod1$residuals
qqPlot(residuos, distribution = "norm", mean = mean(residuos),
sd = sd(residuos))
## [1] 93 92
Aquí podemos ver la distribución de los errores. Básicamente vemos una recta por encima (o debajo) de la cual deberían colocarse los errores, que obviamente, no pasa.
Por las tendencias observadas, es posible pensar que se trata de un modelo exponencial, podemos comprobarlo.
El modelo exponencial se rige por la siguiente expresión:
\[ y = e^{ax+b} = me^{ax} \] Para poder implementar este modelo en R, solo tenemos que aplicar la transformación a los datos, teniendo en cuenta que no existe el logaritmo natural de 0.
mod2 <- lm(log(Casos_Confirmados) ~ Dias,
data = datos_españa[datos_españa$Casos_Confirmados > 0, ])
summary(mod2)
##
## Call:
## lm(formula = log(Casos_Confirmados) ~ Dias, data = datos_españa[datos_españa$Casos_Confirmados >
## 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6795 -1.0578 0.1752 1.1060 1.7823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.008846 0.338741 -8.882 1.36e-13 ***
## Dias 0.193379 0.006012 32.167 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.312 on 81 degrees of freedom
## Multiple R-squared: 0.9274, Adjusted R-squared: 0.9265
## F-statistic: 1035 on 1 and 81 DF, p-value: < 2.2e-16
De este modelo obtenemos p-valores suficientemente baajos como para ver la relevancia de estos parámetros en el modelo. Además, el valor de \(R^2\) es definitivamente mejor que en el caso del modelo lineal. Pero como antes, lo suyo es escribir la ecuación del modelo y luego proceder a representarlo todo.
\[Casos\ Confirmados = 0.0493486 \cdot e^{0.1933787\cdot x}\]
plot(datos_españa$Dias, datos_españa$Casos_Confirmados)
lines(exp(mod2$coefficients[1])*exp(mod2$coefficients[2]*datos_españa$Dias),
col="blue")
Bueno, algo sí que es cierto que ha mejorado, pero ya podemos ver que este modelo seguramente no sea el que buscamos. Podemos hacer un análisis de los errores también. Veamos:
plot(mod2$residuals ~ mod2$fitted.values,
xlab = "Valores ajustados",
ylab = "Residuos del modelo")
residuos = mod2$residuals
qqPlot(residuos, distribution = "norm", mean = mean(residuos),
sd = sd(residuos))
## 34 93
## 24 83
Donde volvemos a ver ciertos patrones, aunque en el análisis de residuos, estos se ajustan mejor en relación al caso anterior. Pero es importante decir que este modelo es simplemente menos pésimo que el lineal, pero pésimo de todos modos. Podemos pensar entonces que se trata de un modelo potencial, vamos a comprobarlo.
En este modelo, aplicaremos el logaritmo a ambas variables, de manera que tendremos algo como:
\[ log(y) = alog(x) + b \Rightarrow y=e^{a\cdot log(x)+b} = m\cdot x^a\] Implementar este modelo en R también es sencillo, y se hace como sigue:
mod3 = lm(log(Casos_Confirmados) ~ log(Dias),
data = datos_españa[datos_españa$Casos_Confirmados > 0, ])
summary(mod3)
##
## Call:
## lm(formula = log(Casos_Confirmados) ~ log(Dias), data = datos_españa[datos_españa$Casos_Confirmados >
## 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8919 -0.9475 0.3698 1.0011 4.7138
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22.6476 1.1733 -19.30 <2e-16 ***
## log(Dias) 7.7885 0.3062 25.44 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.624 on 81 degrees of freedom
## Multiple R-squared: 0.8888, Adjusted R-squared: 0.8874
## F-statistic: 647.1 on 1 and 81 DF, p-value: < 2.2e-16
Como vemos, la \(R^2\) en este caso es un poquito peor que el modelo anterior. Pero lo suyo es representar este modelo de nuevo para verlo.
\[ Casos\ Confirmados = 1.4597421\times 10^{-10}\cdot Dias^{0.1933787} \]
plot(datos_españa$Dias, datos_españa$Casos_Confirmados)
lines(exp(mod3$coefficients[1])*datos_españa$Dias^mod3$coefficients[2],
col = "green")
plot(mod3$residuals ~ mod3$fitted.values,
xlab = "Valores ajustados",
ylab = "Residuos del modelo")
residuos = mod3$residuals
qqPlot(residuos, distribution = "norm", mean = mean(residuos),
sd = sd(residuos))
## 11 12
## 1 2
Si tuviéramos que escoger entre lo que tenemos, me temo que nos tendríamos que quedar con el modelo exponencial. Obviamente esto no es así, ya que existen modelos diferenciales especialmente creados para tratar el estudio pandémico. Este tipo de modelos los veremos más adelante. Pero sí que cabe destacar que podemos complicar el modelo todo lo que queramos. Por ejemplo, podemos inventar uno.
my_model <- lm(log(Casos_Confirmados) ~ Dias + I(Dias^2) + I(Dias^3) + I(Dias^5) + I(Dias^7) + sqrt(Dias), data = datos_españa[datos_españa$Casos_Confirmados > 0, ])
summary(my_model)
##
## Call:
## lm(formula = log(Casos_Confirmados) ~ Dias + I(Dias^2) + I(Dias^3) +
## I(Dias^5) + I(Dias^7) + sqrt(Dias), data = datos_españa[datos_españa$Casos_Confirmados >
## 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23598 -0.07889 0.00744 0.11500 0.57400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.873e+01 1.287e+01 -7.670 4.71e-11 ***
## Dias -1.391e+01 1.675e+00 -8.303 2.89e-12 ***
## I(Dias^2) 2.048e-01 2.375e-02 8.624 7.01e-13 ***
## I(Dias^3) -1.757e-03 2.152e-04 -8.164 5.35e-12 ***
## I(Dias^5) 7.225e-08 1.111e-08 6.505 7.39e-09 ***
## I(Dias^7) -2.049e-12 4.054e-13 -5.053 2.92e-06 ***
## sqrt(Dias) 6.909e+01 8.628e+00 8.008 1.06e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3134 on 76 degrees of freedom
## Multiple R-squared: 0.9961, Adjusted R-squared: 0.9958
## F-statistic: 3247 on 6 and 76 DF, p-value: < 2.2e-16
start_date = ymd("2020-01-22")
end_date = ymd("2020-04-30")
dates = seq(start_date+1, end_date, by = "1 day")
days_since_start = as.numeric(dates - start_date)
new_data = data.frame(Dias = days_since_start)
pred1 = predict(mod1, newdata = new_data)
pred2 = exp(predict(mod2, newdata = new_data))
pred3 = exp(predict(mod3, newdata = new_data))
pred4 = exp(predict(my_model, newdata = new_data))
datos_por_fecha_ts = xts(x = data.frame(Real = c(datos_españa$Casos_Confirmados,rep(NA, length(pred1) - length(datos_españa$Casos_Confirmados))),
Mod_Lin = pred1,
#Mod_Exp = pred2,
Mod_Pot = pred3,
Mod_Mix = pred4),
order.by = dates)
dygraph(datos_por_fecha_ts)
Parece que nuestro modelo mixto, el inventado es el que mejor se ajusta a la pandemia, aunque obviamente aún comete cierto error. A partir de los datos reales es normal que ya se descuadre, porque está ajustado según estos datos. Más adelante utilizaremos los modelos que se suelen utilizar en el estudio de epidemias y podremos comparar.
Para poder realizar un estudio de la homogeneidad de la infección, vamos a cargar de nuevo el csv, para que no interfiera con nuestros análisis previos. Por ello, le asignaremos otro nombre.
covid19 = read.csv("covid_19_clean_complete.csv")
Para calcular el número de habitantes por país, vamos a utilizar otra librería muy útil que es wbstats (ponemos entre 2018 y 2019 porque aún no dispone de los datos para 2020).
pop_data <- wb(indicator = "SP.POP.TOTL", startdate = 2018, enddate = 2019)
pop_data %>% head() %>% kable()
| iso3c | date | value | indicatorID | indicator | iso2c | country | |
|---|---|---|---|---|---|---|---|
| 2 | ARB | 2018 | 419790588 | SP.POP.TOTL | Population, total | 1A | Arab World |
| 4 | CSS | 2018 | 7358965 | SP.POP.TOTL | Population, total | S3 | Caribbean small states |
| 6 | CEB | 2018 | 102511922 | SP.POP.TOTL | Population, total | B8 | Central Europe and the Baltics |
| 8 | EAR | 2018 | 3249140605 | SP.POP.TOTL | Population, total | V2 | Early-demographic dividend |
| 10 | EAS | 2018 | 2328220870 | SP.POP.TOTL | Population, total | Z4 | East Asia & Pacific |
| 12 | EAP | 2018 | 2081651801 | SP.POP.TOTL | Population, total | 4E | East Asia & Pacific (excluding high income) |
Con esto, podemos calcular la población de distintos países, por ejemplo:
pop_data[pop_data$country == "Spain", ]
## iso3c date value indicatorID indicator iso2c country
## 454 ESP 2018 46796540 SP.POP.TOTL Population, total ES Spain
pop_data[pop_data$country == "Italy", ]
## iso3c date value indicatorID indicator iso2c country
## 288 ITA 2018 60421760 SP.POP.TOTL Population, total IT Italy
pop_data[pop_data$country == "France", ]
## iso3c date value indicatorID indicator iso2c country
## 232 FRA 2018 66977107 SP.POP.TOTL Population, total FR France
Obviamente, habrá países de los que no dispongamos de los datos, es por esto que podemos hacer un pequeño data wrangling para eliminarlos:
paises = unique(covid19$Country.Region)
covid19.limpio = c()
for (i in 1:length(paises)){
if(length(which(paises[i] %in% pop_data$country)) > 0){
covid19.limpio = rbind(covid19.limpio, covid19[covid19$Country.Region==paises[i],])
}
}
Podemos realizar una pequeña representación de los datos que tenemos hasta ahora. Para ello podemos utilizar la función aggregate.
covid19.limpio$Date=as.Date(as.character(covid19.limpio$Date),"%m/%d/%Y")
infectados.totales.por.dia = aggregate(covid19.limpio$Confirmed ~
covid19.limpio$Date, FUN = sum)
fallecidos.totales.por.dia = aggregate(covid19.limpio$Deaths ~
covid19.limpio$Date, FUN = sum)
recuperados.totales.por.dia = aggregate(covid19.limpio$Recovered ~
covid19.limpio$Date, FUN = sum)
tabla.totales = data.frame(infectados.totales.por.dia[, 1],
infectados.totales.por.dia[, 2],
fallecidos.totales.por.dia[, 2],
recuperados.totales.por.dia[, 2])
names(tabla.totales) = c("Fecha", "Infectados", "Fallecidos", "Recuperados")
fig <- plot_ly(
tabla.totales,
x = ~tabla.totales$Fecha,
y = ~tabla.totales$Fallecidos,
name = 'Fallecidos',
type = 'scatter',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~tabla.totales$Infectados,
name = 'Infectados',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~tabla.totales$Recuperados,
name = 'Recuperados',
mode = 'lines')
#Damos formato a la figura
fig <- fig %>% layout(title = 'Análisis temporal' ,yaxis = list(title = 'Número de Casos'), xaxis = list(title = 'Fecha'))
fig
Para poder hacer el test de \(\chi^2\) en R, primero debemos fijar una fecha, para luego realizar los cálculos pertinentes.
fecha = "0020-03-30"
confirmados.por.pais = aggregate(covid19.limpio$Confirmed[covid19.limpio$Date == fecha] ~ covid19.limpio$Country.Region[covid19.limpio$Date == fecha], FUN = sum)
names(confirmados.por.pais)=c("Pais","Confirmados")
paises = unique(covid19.limpio$Country.Region)
suma.total.habitantes = sum(pop_data[pop_data$country %in% paises,]$value)
numero.total.infectados = sum(confirmados.por.pais$Confirmados)
Para poder realizar creamos una tabla nueva donde aparezcan las frecuencias estimadas.
tabla.infectados.paises = c()
for (i in 1:length(paises)){
habitantes = pop_data[pop_data$country == paises[i],]$value
confirmados = confirmados.por.pais$Confirmados[confirmados.por.pais$Pais==paises[i]]
confirmados.estimados = numero.total.infectados*habitantes/suma.total.habitantes
tabla.infectados.paises = rbind(tabla.infectados.paises,
c(confirmados, confirmados.estimados))
}
tabla.infectados.paises = as.data.frame(tabla.infectados.paises)
tabla.infectados.paises = data.frame(paises, tabla.infectados.paises)
names(tabla.infectados.paises) = c("pais", "infectados", "infectados.estimados")
chisq.test(tabla.infectados.paises$infectados,
p = tabla.infectados.paises$infectados.estimados/sum(tabla.infectados.paises$infectados))
## Warning in chisq.test(tabla.infectados.paises$infectados,
## p = tabla.infectados.paises$infectados.estimados/
## sum(tabla.infectados.paises$infectados)): Chi-squared approximation may be
## incorrect
##
## Chi-squared test for given probabilities
##
## data: tabla.infectados.paises$infectados
## X-squared = 5625480, df = 157, p-value < 2.2e-16
Vemos que nos salen algunos problemas, para solucionarlos hagamos lo siguiente:
paises.con.problemas = which(tabla.infectados.paises$infectados.estimados < 5)
paises[paises.con.problemas]
## [1] Liechtenstein Monaco San Marino
## 185 Levels: Afghanistan Albania Algeria Andorra Angola ... Zimbabwe
tabla.infectados.paises2 = tabla.infectados.paises[-paises.con.problemas,]
pais.añadir = data.frame("problemas", sum(tabla.infectados.paises[
tabla.infectados.paises$pais %in% paises[paises.con.problemas],]
$infectados), sum(tabla.infectados.paises[tabla.infectados.paises$pais %in% paises[paises.con.problemas],]$infectados.estimados))
names(pais.añadir) = names(tabla.infectados.paises2)
pais.añadir
## pais infectados infectados.estimados
## 1 problemas 341 9.464223
tabla.infectados.paises2 = rbind(tabla.infectados.paises2, pais.añadir)
chisq.test(tabla.infectados.paises2$infectados,
p = tabla.infectados.paises2$infectados.estimados/sum(tabla.infectados.paises2$infectados))
##
## Chi-squared test for given probabilities
##
## data: tabla.infectados.paises2$infectados
## X-squared = 5617599, df = 155, p-value < 2.2e-16
Al obtener un p-valor tan pequeño, podemos concluir que el virus no se expandió del mismo modo por todos los países el día 30 de marzo, que era el resultado esperado.
datos_CCAA = read.csv("DatosCCAA.csv")
datos_CCAA %>% head() %>% kable()
| fecha | cod_ine | CCAA | total | tipo | Pob_CCAA_2019 | UCIS_Públicos | UCIS_Privados | UCI_Total |
|---|---|---|---|---|---|---|---|---|
| 2020-02-27 | 1 | Andalucía | 1 | casos | 8414240 | 572 | 162 | 734 |
| 2020-02-28 | 1 | Andalucía | 6 | casos | 8414240 | 572 | 162 | 734 |
| 2020-02-29 | 1 | Andalucía | 8 | casos | 8414240 | 572 | 162 | 734 |
| 2020-03-01 | 1 | Andalucía | 12 | casos | 8414240 | 572 | 162 | 734 |
| 2020-03-02 | 1 | Andalucía | 12 | casos | 8414240 | 572 | 162 | 734 |
| 2020-03-03 | 1 | Andalucía | 13 | casos | 8414240 | 572 | 162 | 734 |
A priori estos datos son complicados de entender. Por ello vamos a hacer alguna transformación.
datos_CCAA %<>%
select(fecha, tipo, CCAA, total)
datos_CCAA %<>%
pivot_wider(names_from = CCAA, values_from = total)
datos_CCAA %>% head() %>% kable()
| fecha | tipo | Andalucía | Aragón | Asturias | Baleares | Canarias | Cantabria | Castilla-La Mancha | Castilla y León | Cataluña | C. Valenciana | Extremadura | Galicia | Madrid | Murcia | Navarra | País Vasco |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2020-02-27 | casos | 1 | 0 | 0 | 1 | 6 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 4 | 0 | 0 | 0 |
| 2020-02-28 | casos | 6 | 1 | 0 | 1 | 6 | 0 | 0 | 2 | 3 | 8 | 0 | 0 | 5 | 0 | 0 | 0 |
| 2020-02-29 | casos | 8 | 1 | 0 | 2 | 6 | 0 | 0 | 2 | 5 | 10 | 0 | 0 | 8 | 0 | 0 | 2 |
| 2020-03-01 | casos | 12 | 0 | 1 | 2 | 7 | 1 | 1 | 3 | 6 | 15 | 4 | 0 | 10 | 0 | 1 | 3 |
| 2020-03-02 | casos | 12 | 0 | 1 | 2 | 7 | 10 | 3 | 3 | 15 | 15 | 6 | 0 | 29 | 0 | 1 | 9 |
| 2020-03-03 | casos | 13 | 0 | 1 | 2 | 7 | 10 | 7 | 8 | 15 | 15 | 6 | 0 | 49 | 0 | 1 | 13 |
#str(datos_CCAA)
Ahora para poder representar esto, hagamos lo siguiente:
casos_ccaa = datos_CCAA %>%
filter(datos_CCAA$tipo == "casos")
casos_ccaa %>% head() %>% kable()
| fecha | tipo | Andalucía | Aragón | Asturias | Baleares | Canarias | Cantabria | Castilla-La Mancha | Castilla y León | Cataluña | C. Valenciana | Extremadura | Galicia | Madrid | Murcia | Navarra | País Vasco |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2020-02-27 | casos | 1 | 0 | 0 | 1 | 6 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 4 | 0 | 0 | 0 |
| 2020-02-28 | casos | 6 | 1 | 0 | 1 | 6 | 0 | 0 | 2 | 3 | 8 | 0 | 0 | 5 | 0 | 0 | 0 |
| 2020-02-29 | casos | 8 | 1 | 0 | 2 | 6 | 0 | 0 | 2 | 5 | 10 | 0 | 0 | 8 | 0 | 0 | 2 |
| 2020-03-01 | casos | 12 | 0 | 1 | 2 | 7 | 1 | 1 | 3 | 6 | 15 | 4 | 0 | 10 | 0 | 1 | 3 |
| 2020-03-02 | casos | 12 | 0 | 1 | 2 | 7 | 10 | 3 | 3 | 15 | 15 | 6 | 0 | 29 | 0 | 1 | 9 |
| 2020-03-03 | casos | 13 | 0 | 1 | 2 | 7 | 10 | 7 | 8 | 15 | 15 | 6 | 0 | 49 | 0 | 1 | 13 |
fig <- plot_ly(
casos_ccaa,
x = ~casos_ccaa$fecha,
y = ~casos_ccaa$Andalucía,
name = 'Andalucía',
type = 'scatter',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Aragón,
name = 'Aragón',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Asturias,
name = 'Asturias',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Baleares,
name = 'Baleares',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Canarias,
name = 'Canarias',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Cantabria,
name = 'Canatabria',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$`Castilla-La Mancha`,
name = 'Castilla La Mancha',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$`Castilla y León`,
name = 'Castilla y León',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Cataluña,
name = 'Cataluña',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$`C. Valenciana`,
name = 'C. Valenciana',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Extremadura,
name = 'Extremadura',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Galicia,
name = 'Galicia',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Madrid,
name = 'Madrid',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Murcia,
name = 'Murcia',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$Navarra,
name = 'Navarra',
mode = 'lines')
fig <- fig %>% add_trace(
y = ~casos_ccaa$`País Vasco`,
name = 'País Vasco',
mode = 'lines')
#Damos formato a la figura
fig <- fig %>% layout(title = 'Análisis por CCAA' ,yaxis = list(title = 'Número de Casos'), xaxis = list(title = 'Fecha'))
fig
Para poder hacer el test de \(\chi^2\), debemos fijar de nuevo una fecha, para ello hacemos:
datos_ccaa2 = read.csv("DatosCCAA.csv")
fecha = "2020-03-15"
confirmados.por.ccaa = aggregate(datos_ccaa2$total[datos_ccaa2$fecha == fecha] ~ datos_ccaa2$CCAA[datos_ccaa2$fecha == fecha], FUN = sum)
names(confirmados.por.ccaa)=c("CCAA","Confirmados")
confirmados.por.ccaa %>% head() %>% kable()
| CCAA | Confirmados |
|---|---|
| Andalucía | 443 |
| Aragón | 154 |
| Asturias | 138 |
| Baleares | 29 |
| C. Valenciana | 414 |
| Canarias | 110 |
poblacion <- read.csv("Fichero_poblaciónCCAA.csv")
ccaa = unique(datos_ccaa2$Pob_CCAA_2019)
ccaa2 = unique(poblacion$CCAA)
ccaa3 = unique(datos_ccaa2$CCAA)
suma.total.habitantes = as.numeric(sum(ccaa))
numero.total.infectados = sum(confirmados.por.ccaa$Confirmados)
tabla.infectados.ccaa = c()
for (i in 1:length(ccaa2)){
habitantes = as.numeric(poblacion[poblacion$CCAA == ccaa2[i],]$Pob_CCAA_2019)
confirmados = confirmados.por.ccaa$Confirmados[confirmados.por.ccaa$CCAA==ccaa3[i]]
confirmados.estimados = numero.total.infectados*habitantes/suma.total.habitantes
tabla.infectados.ccaa = rbind(tabla.infectados.ccaa,
c(confirmados, confirmados.estimados))
}
tabla.infectados.ccaa = as.data.frame(tabla.infectados.ccaa)
tabla.infectados.ccaa = data.frame(ccaa2, tabla.infectados.ccaa)
names(tabla.infectados.ccaa) = c("pais", "infectados","infectados.estimados")
tabla.infectados.ccaa %>% head() %>% kable()
| pais | infectados | infectados.estimados |
|---|---|---|
| ANDALUCÍA | 443 | 1397.42698 |
| ARAGÓN | 154 | 219.10628 |
| ASTURIAS | 138 | 169.86541 |
| ILLES BALEARS | 29 | 190.90095 |
| CANARIAS | 110 | 357.63229 |
| CANTABRIA | 51 | 96.50474 |
chisq.test(tabla.infectados.ccaa$infectados,
p = tabla.infectados.ccaa$infectados.estimados/sum(tabla.infectados.ccaa$infectados))
##
## Chi-squared test for given probabilities
##
## data: tabla.infectados.ccaa$infectados
## X-squared = 8386.1, df = 15, p-value < 2.2e-16
Por tanto, podemos concluir, como era esperado, que la pandemia no se propagó de igual manera para todas las CCAA a día 15 de marzo de 2020.
Como sabemos, el modelo SIR, define las siguientes variables:
En el modelo, se definen además los siguientes parámetros:
De esta manera, tendremos una tasa de transición, definida por:
\[ \beta \cdot \frac{I(t)}{N} \] Y una tasa de transición, definida directamente por \(\gamma\). De modo que se definen las ecuaciones diferenciales del modelo SIR, que son las que trataré de implementar, de la siguiente forma:
\[ \frac{dS}{dt}=-\frac{\beta I}{N}\cdot S \], \[ \frac{dI}{dt} = \frac{\beta I}{N}\cdot S - \gamma I\], \[\frac{dR}{dt}= \gamma I\]
La resolución de estas ecuaciones no es el objetivo de este estudio, pero se puede investigar en la literatura. Para empeazar, cargamos los datos (que ya hicimos anteriormente), pero esta vez usaremos date. Por evitar posibles problemas de interacción, volveré a cargar los datos, aunque no es necesario.
covid_19 = read.csv("covid_19_clean_complete.csv")
covid_19$Date = as.Date(as.character(covid_19$Date), "%m/%d/%Y")
pop_data <- wb(indicator = "SP.POP.TOTL", startdate = 2018, enddate = 2019)
Nos interesan sobretodo los datos referidos a España, para ello:
covid19_España = covid_19[covid_19$Country.Region == 'Spain', ]
infectados_por_dia = aggregate(covid19_España$Confirmed ~ covid19_España$Date,
FUN = sum)
fallecidos_por_dia = aggregate(covid19_España$Deaths ~ covid19_España$Date,
FUN = sum)
infectados_por_dia2 = infectados_por_dia[, 2] + fallecidos_por_dia[, 2]
#Sumamos infectados y fallecidos porque el modelo SIR no considera fallecidos.
recuperados_por_dia = aggregate(covid19_España$Recovered ~covid19_España$Date,
FUN = sum)
Los susceptibles se calculan del siguiente modo:
habitantes = pop_data$value[pop_data$country=="Spain"]
susceptibles.por.dia = habitantes - infectados_por_dia2 - recuperados_por_dia[,2]
Por último, creamos una tabla con toda la información:
tabla.España = data.frame(unique(covid19_España$Date), susceptibles.por.dia,
infectados_por_dia2, recuperados_por_dia[,2])
names(tabla.España) = c("Fecha", "Susceptibles", "Infectados", "Recuperados")
Para estimar el valor de R0, hagamos lo siguiente:
x = tabla.España$Recuperados
y = habitantes * log(tabla.España$Susceptibles)
summary(lm(y ~ x))
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40930 -2940 8669 8678 26980
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.265e+08 1.983e+03 416786.51 <2e-16 ***
## x -3.853e+00 6.367e-02 -60.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16520 on 91 degrees of freedom
## Multiple R-squared: 0.9757, Adjusted R-squared: 0.9755
## F-statistic: 3661 on 1 and 91 DF, p-value: < 2.2e-16
El valor de R0 es, por tanto, 3.853, con un \(R^2\) de 0.97, lo cual parece bueno. Así pues, la estimación de R0 es:
(estimacion.R0 = -summary(lm(y ~ x))$coefficients[2])
## [1] 3.852802
dia.ultimo = length(tabla.España[, 1])
exp(estimacion.R0 * tabla.España$Recuperados[dia.ultimo]/habitantes)
## [1] 1.007375
Concluimos pues, que como el valor de R0 supera al valor del cálculo exponencial, el coronavirus, hasta la fecha de los datos en la que se están utilizando para este análisis, se seguirá propagando.
Una serie temporal puede considerarse un proceso estocástico en el que se estudian la evolución en el tiempo de una sucesión de variables aleatorias. Dichas variables aleatorias normalmente están correlacionadas o existe una dependencia lineal entre las mismas.
En una serie temporal Y(t), se intenta identificar 4 componentes:
Hay varios modelos de series temporales, como el aditivo o el multiplicativo. De esto se puede investigar en la literatura. Normalmente, nos tendremos que fijar en la amplitud de la componente estacional, que será la que determinará de qué tipo es la serie temporal de la que estamos tratando.
Hacer predicciones con este tipo de modelos no es siempre fácil, si no que se utilizan varios modelos, por ejemplo el modelo de Holt-Winters, que será el que utilizaremos.
Como siempre, tenemos que empezar por la carga de datos, y el filtrado para quedarnos únicamente con los datos de España:
covid_19 = read.csv("covid_19_clean_complete.csv")
covid_19$Date = as.Date(as.character(covid_19$Date), "%m/%d/%Y")
covid19_España = covid_19[covid_19$Country.Region == 'Spain', ]
infectados_por_dia = aggregate(covid19_España$Confirmed ~ covid19_España$Date,
FUN = sum)
fallecidos_por_dia = aggregate(covid19_España$Deaths ~ covid19_España$Date,
FUN = sum)
recuperados_por_dia = aggregate(covid19_España$Recovered ~covid19_España$Date,
FUN = sum)
Creamos una tabla con estos datos:
tabla.España = data.frame(unique(covid19_España$Date), infectados_por_dia[,2],
fallecidos_por_dia[,2], recuperados_por_dia[,2])
names(tabla.España) = c("Fecha", "Infectados", "Fallecidos", "Recuperados")
Ahora que tenemos los datos preparados, podemos hacer un estudio para los infectados, por ejemplo. Para poder hacerlo, necesitamos modelizar los infectados como serie temporal, esto es:
infectados = ts(tabla.España$Infectados, frequency = 7, start = c(1,3))
#infectados
Donde hemos definido una estacionalidad semanal (frequency = 7), y que el primer día fue miércoles, que es el primer dato que tenemos. Usando la librería forecast, podemos hacer representaciones como la siguiente:
autoplot(infectados)
Viendo esto, supondremos en primer lugar que nuestra serie es aditiva, por ello, hagamos lo siguiente:
components = decompose(infectados, type = "additive")
#components$seasonal
autoplot(components)
Aquí tenemos las 4 componentes de las que hablábamos, donde podemos ver en la aleatoria, que en torno a la semana número 10 empieza a crecer. Esto es malo, pues lo que realmente nos interesa es que esta componente se mantenga constante a lo largo del tiempo. Para poder ajustar nuestra serie temporal, debemos eliminar la componente estacional, por tanto, hagamos lo siguiente:
infectados_ajustado = infectados - components$seasonal
autoplot(infectados_ajustado)
Con este modelo ajustado, podemos hacer nuestras predicciones. En este caso, como ya dije, utilizaré el suavizado de Holt-Winters.
(prediccion_infectados = HoltWinters(infectados_ajustado, gamma = FALSE))
## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = infectados_ajustado, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.7335536
## beta : 0.769321
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 212663.450
## b 4096.022
El valor de \(\alpha\) es alto, lo que nos dice que las predicciones usan el nivel actual de la serie de infecyados con un peso elevado. El valor de \(\beta\) nos habla de la pendiente de la tendencia de la predicción en un día particular. Podemos ver el error cometido en nuestra predicción de la siguiente manera.
prediccion_infectados$SSE
## [1] 108479941
sqrt(prediccion_infectados$SSE)
## [1] 10415.37
plot(prediccion_infectados)
Pero lo interesante es hacer predicciones más allá del último día del que disponemos. Para ello utilizaremos la función forecast.HoltWinters. Aquí, el parámetro h indicará hasta dónde queremos realizar la predicción (pondremos, por ejemplo, 7 semanas).
(prediccion_infectados_semanal = forecast:::forecast.HoltWinters(prediccion_infectados, h=7))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 14.57143 216759.5 215356.9 218162.0 214614.5 218904.5
## 14.71429 220855.5 218557.5 223153.5 217341.0 224370.0
## 14.85714 224951.5 221472.6 228430.4 219631.0 230272.0
## 15.00000 229047.5 224180.7 233914.3 221604.4 236490.7
## 15.14286 233143.6 226718.4 239568.7 223317.1 242970.0
## 15.28571 237239.6 229106.5 245372.7 224801.1 249678.1
## 15.42857 241335.6 231359.1 251312.2 226077.8 256593.4
Aquí podemos ver los intervalos de confianza. Para poder representarlo, hacemos lo mismo y así ver las bandas de los intervalos de confianza.
autoplot(prediccion_infectados_semanal)
Para comprobar el modelo, usaremos el test de Ljung-Box y ver si las autocorrelaciones de los errores de la predicción de nuestra serie son diferentes de 0 o no. Para que el modelo sea válido, deben ser cero. El correlograma se define del siguiente modo:
ggAcf(prediccion_infectados_semanal$residuals[3:92], lag.max = 7)
Box.test(prediccion_infectados_semanal$residuals, lag = 7, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: prediccion_infectados_semanal$residuals
## X-squared = 5.4665, df = 7, p-value = 0.6032
También podemos comprobar la normalidad de los residuos con el test de Saphiro-Wilks:
shapiro.test(prediccion_infectados_semanal$residuals)
##
## Shapiro-Wilk normality test
##
## data: prediccion_infectados_semanal$residuals
## W = 0.86898, p-value = 1.868e-07
Vemos, con este p-valor tan pequeño, que no siguen una distribución normal y por tanto, el modelo no sería adecuado. Podemos hacer esto también para los fallecidos.
fallecidos = ts(tabla.España$Fallecidos, frequency = 7, start = c(1,3))
autoplot(fallecidos)
components = decompose(fallecidos, type = "additive")
#components$seasonal
autoplot(components)
fallecidos_ajustado = fallecidos - components$seasonal
autoplot(fallecidos_ajustado)
(prediccion_fallecidos = HoltWinters(fallecidos_ajustado, gamma = FALSE))
## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = fallecidos_ajustado, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.8473805
## beta : 0.7639737
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 22141.3889
## b 419.7952
sqrt(prediccion_fallecidos$SSE)
## [1] 887.4306
plot(prediccion_fallecidos)
(prediccion_fallecidos_semanal = forecast:::forecast.HoltWinters(prediccion_fallecidos, h = 7))
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 14.57143 22561.18 22441.69 22680.68 22378.43 22743.94
## 14.71429 22980.98 22766.08 23195.88 22652.31 23309.65
## 14.85714 23400.77 23066.55 23735.00 22889.62 23911.93
## 15.00000 23820.57 23348.53 24292.61 23098.65 24542.49
## 15.14286 24240.36 23614.67 24866.06 23283.45 25197.28
## 15.28571 24660.16 23866.63 25453.69 23446.57 25873.75
## 15.42857 25079.96 24105.60 26054.31 23589.80 26570.11
autoplot(prediccion_fallecidos_semanal)
ggAcf(prediccion_fallecidos_semanal$residuals[3:92], lag.max=7)
Box.test(prediccion_fallecidos_semanal$residuals, lag = 7, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: prediccion_fallecidos_semanal$residuals
## X-squared = 16.249, df = 7, p-value = 0.02294
shapiro.test(prediccion_fallecidos_semanal$residuals)
##
## Shapiro-Wilk normality test
##
## data: prediccion_fallecidos_semanal$residuals
## W = 0.80548, p-value = 1.308e-09
Podemos concluir por tanto, que el modelo no es el adecuado. Quizá esto sea porque el modelo es multiplicativo u otro tipo más complejo.
Tras el análisis realizado, podemos observar que el avance de la pandemia se ha ido frenando, en parte gracias a las medidas de confinamiento que se han tomado en la mayor parte de países afectados. A pesar de que las predicciones realizadas en este análisis son básicas, nos han servido para tener una idea general de lo que puede llegar a pasar.
Los datos de este estudio se pueden ir actualizando según la base de datos dada en el enlace mencionado al principio del análisis.